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The thermal behavior of free and alumina-supported iron-carbon nanoparticles is investigated 
via molecular dynamics simulations, in which the effect of the substrate is treated with a simple 
Morse potential fitted to ab initio data. We observe that the presence of the substrate raises the 
melting temperature of medium and large Fe\- x C x nanoparticles (x = — 0.16, N = 80 — 1000, non- 
magic numbers) by 40-60 K; it also plays an important role in defining the ground state of smaller 
Fe nanoparticles (N = 50 — 80). The main focus of our study is the investigation of Fe-C phase 
diagrams as a function of the nanoparticle size. We find that as the cluster size decreases in the 
1.1-1.6-nm-diameter range the eutectic point shifts significantly not only toward lower temperatures, 
as expected from the Gibbs- Thomson law, but also toward lower concentrations of C. The strong 
dependence of the maximum C solubility on the Fe-C cluster size may have important implications 
for the catalytic growth of carbon nanotubes by chemical vapor deposition. 
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I. INTRODUCTION 

Catalytic chemical vapor deposition (CVD) is a widely 
used method for the production of carbon nanotubes 
(CNTs) by decomposition of hydrocarbons (such as CH4, 
C2H2 etc.) or carbon monoxide on snpported metal cat- 
alysts (Fe, Ni, Co, FeMo, etc.) 0, Still i]. Despite 
numerous studies, the growth mechanism of nanotubes is 
still not well understood. Among the most studied fac- 
tors that control the growth process are the kinetics of 
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The vapor-liquid-solid (VLS) model [ljj for the CNT 
growth by CVD implies that the catalytic particle should 
be in a liquid state which allows rapid diffusion of car- 
bon atoms throughout the particle. The bulk diffusion of 
carbon through the metal nanoparticle is driven by con- 
centration gradients [2(j, and it is considered to be the 
rate-limiting step in the growth of filaments or carbona- 
ceous deposits [UHII- The activation energies (~1.2-1.8 
eV) measured for thermal CVD growth of nanofibers or 
nanotubes are consistent with those for the carbon dif- 
fusion through the corresponding metals 0, E 0, I22I l23l|. 
hence further supporting the bulk diffusion VLS model. 
Another mechanism, the surface-mediated carbon trans- 
port model, has been proposed E 0|- The low tem- 
perature nanotube synthesis by plasma-enhanced CVD 
E E El El HI El implies that the catalyst could be 
in a solid state. However, the temperature of the ac- 



tive nanoparticles is extremely difficult to measure dur- 
ing the growth process. In fact, the bombardment of 
energetic species in plasma and the exothermic dehydro- 
genation reactions of hydrocarbon can increase the local 
temperature of the catalyst, even if the substrate is kept 
in thermal equilibrium. This phenomenon can promote 
the surface melting of the nanoparticle and, concurrently, 
facilitate the growth of the nanotube [H EH • 

Because the overall catalytic capability of nanoparti- 
cles strongly depends on whether they are in the liquid 
or solid states, their thermal behavior has been exten- 
sively investigated with experimental I 111 . Il2t I26L mm 
29|, |M Jll|jM_ai 



and theoretical means [13|, [14|, [33|, [34|, 
[40| . Under CVD experimental condi- 
tions, the melting temperatures of catalytic particles are 
strongly reduced because of the dissolved carbon (liq- 
uidus and solidus slopes in the metal-C phase diagram) 
and the relatively high surface energy with respect to the 
bulk materials (Gibbs-Thomson phenomenon). Unusu- 
ally low melting points of 600-700°C have been observed 
for oversaturated solutions of carbon (up to 50 at.% at 
700°C) in Fe, Ni, and Co metals [3l|. In these cases, 
the fluidization of the metal catalytic particles at low 
or moderate temperature was attributed to the creation 
of highly dispersed unstable solutions oversaturated with 
carbon, with concentrations well above the limit of the 
stable carbide [32|, |4l| . 

With so many factors influencing the catalytic growth 
of CNT, the role of the substrate on the thermal proper- 
ties of the particle is often overlooked. In fact, there are 
experimental indications that the presence of substrate 
could be an important factor in the thermodynamics of 
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the particle and in regulating the growth of nanotubes 
0, 113, Gil- Most of the theoretical research has been 
focused on the melting behavior of free-standing pure 
[lllMIMIl or bimetallic [MS 11 IIS III clusters. A 
phase diagram of free-standing Fe-C clusters of fixed size 
(~ 2.4 nm) has been recently calculated [HI, |4^|. There 
have been investigations performed on supported cata- 
lysts [H, [H, |H, |3] , showing that the cluster-substrate 
interaction strongly affects the melting temperature of 
particles and thus could influence the nanotube growth 
rate. A more general picture, capturing the particle- 
substrate interaction and size effects on the phase dia- 
gram, is still lacking and is the main subject of our study. 

In this paper we investigate thermal behaviors for free 
and alumina-supported Fe-C nanoparticles with molec- 
ular dynamics (MD) simulations. In Section II we de- 
scribe the ab initio development of a simple interaction 
between Fe and AI2O3 which is essential for the calcula- 
tions of supported clusters of reasonable sizes. In Section 
III, by using MD simulations we explore the thermody- 
namics of Fe nanoclusters and the phase diagrams of Fe- 
C binary nanoparticles focusing on the effects of cluster's 
size and interfacial interactions on the thermal properties 
(Fe-C up to ~16% carbon concentration). In the same 
Section we show that not only the eutectic temperature 
but also the eutectic composition are size- and substrate- 
dependent. Section IV is devoted to the exploration of 
the peculiar characteristics of very small clusters. Con- 
clusions are described in Section V. 



II. THE FE-AL2O3 POTENTIAL 

In this section we model the many-body Fe-A^Os in- 
teraction with a simple classical potential. If the cal- 
culation of this interaction energy involved the sum- 
mation over all substrate atoms (e.g. using modified 
charge transfer potential [45|, |4(|) the simulation would 
be computationally very expensive (a 200-atom nanopar- 
ticle would require a 1000-atom substrate patch). Ideally, 
a suitable potential for MD simulations would depend on 
three integral variables, or even one (the distance from 
an iron atom to the surface) if the potential corrugation 
is small. The validity of such simplification critically de- 
pends on how weak the Fe-A^Oa interaction is compared 
to those of Fe-Fe and AI2O3-AI2O3. 

There are three possible terminations for (0001) a- 
A1 2 3 surfaces [H IH, H HI : stoichiometric Al- 
terminated, Al-Al-terminated (two top Al layers), and 
O-terminated. These are depicted in Figure QJa) . The- 
oretical calculations have predicted that the most stable 
surface is the stoichiometric one, terminated by a sin- 
gle layer of Al [47| . This type of surface is also believed 
to be the most often observed under ultra high vacuum 
conditions [H [4t| [5(| • Therefore in this work we will re- 
strict our analysis of Fe over AI2O3 using only the stable 
stoichiometric Al-termination. 



A. Details of Calculation 

To develop the interaction between Fe clusters and 
AI2O3 surfaces, we perform density functional theory 
ab initio calculations with Vienna Ab-initio Simulation 
Package VASP [5l|, |52| with projector augmented waves 
(PAW) [H, HiJ and exchange-correlation functionals as 
parameterized by Perdew, Burke, and Ernzerhof (PBE) 
[55j | for the generalized gradient approximation (GGA). 
Simulations are carried out at zero temperature, with 
spin polarization and without zero-point motion. We 
use an energy cutoff of 500 eV and a 6 x 6 x 1 k-point 
Monkhorst-Pack mesh (5(| . The force tolerance for struc- 
tural relaxation is set to 0.1 eV/A. The unit cell is hexag- 
onal with lattice parameters a = 4.767 and c = 29.143 
A. The length of the lattice vector normal to surface is 
kept large enough to minimize the neighboring supercell 
interaction. Vertically, there is at least 12 A of empty 
space. In addition, a dipole layer is applied in z-direction 
to minimize electrostatic effects. A schematic of the unit 
cell is shown in Figure W[b) . 




FIG. 1: (color online), (a) Possible surface terminations of a- 
AI2O3 (side view), (b) Schematic of the hexagonal unit cell, 
(c) The two different high symmetry positions "1" and "2" of 
the inner layer of adsorbed iron (top view), each Fe layer has 
four atoms. 



B. Results and Discussion 

The calculated lattice parameters of the bulk alumina 
(space-group # 167) with a hexagonal unit cell (hRSO), 
a = 4.767 A and c = 12.994 A, are in good agree- 
ment with the experimental values, a exp = 4.742 A and 
c exp = 12.919 A [13, HI- Our a-Al 2 3 slab consists of 
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four oxygen layers as shown in Figure [T^b) . The struc- 
ture is relaxed from the bulk a-Al2 03 configuration while 
keeping one oxygen and one aluminum layers at the bot- 
tom frozen. The surface undergoes considerable relax- 
ations: the interlayer spacings for the top four layers be- 
come d = 0.117,0.899, 1.019 and 0.264 A, differing from 
the corresponding bulk values by -87.3%, 7.8%, -47.1%, 
and 22.1 %, respectively. These values compare well with 
the previous 18 oxygen-layer simulations (- 87.4%, 3.1%, - 
41.7%, 18.9%) [47j, suggesting that an alumina slab with 
four layers of oxygen is thick enough to be a realistic 
model of the AI2O3 surface. On the top of the relaxed a- 
AI2O3 substrate we place a few close-packed layers of iron 
while keeping the in-plane lattice vectors of Fe identical 
to the bulk value of alumina. Because of the fortunate 
match between the natural Fe and the a-A^O^ surface 
lattice spacings (within 3%), the effect of the interface 
strain energy on the total binding is insignificant. Close- 
packed iron layers can be put in different ways on top 
of alumina. Thus, we consider the two high symmetry 
positions labeled as "1" and "2" in Figure [He). 
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FIG. 2: (color online). The definition of the binding and 
deformation energies. 



The binding energy between one Fe column (one Fe 
per layer) and AI2O3 is defined as 



E h = - \E 



Fe-Al 2 3 
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(1) 



where Ep e -Ai 2 o 3 is the total energy of the optimized 
system, E re i axe MAi 2 o 3 ] an( i Ep e are the energies of the 
relaxed AI2O3 slab and Fe layers, respectively (note that 
the factor 1/4 appears because there are four Fe per 
adsorbed layer). The presence of adsorbed Fe strongly 
modifies the alumina surface, and the deformation of the 
surface may affect the total binding considerably. To 
evaluate this effect, we define the substrate deformation 



where Edeformed[Ai 2 o 3 } is the energy of an artificially iso- 
lated alumina slab with the same geometry as that of the 
optimized Fe-Al203 system but without adsorbed iron. 
Figure shows a schematic of the definitions of the two 
energies, Ei, and Edeform- In our calculations only the 
surface is allowed to relax: iron atoms are kept at their 
ideal positions because we assume that a reliable Fe-Fe 
interaction (e.g. Born-Mayer potential [H^, [6(3]) will be 
able to properly describe their energetics. 

To determine how many layers of Fe must be included 
in the calculations, we evaluate the Fe-Al203 interaction 
for systems with different number of Fe layers situated at 
the same distance above the alumina surface. As sum- 
marized in Table HI the binding energy converges reason- 
ably well once the iron film contains three or more layers. 
Therefore, we limit our calculations to the absorption of 
only three Fe layers. In addition, due to the corrugation 
of the surface, the binding energy depends on where the 
adsorbants are positioned in the x-y plane. Figure [ljc) 
shows the two high symmetry configurations of the inner 
layer of adsorbed film. The corrugation, calculated as 
the difference between the binding energies of the most 
and least stable configurations (position "1" and "2" in 
Figure QJc), respectively), is less than 15% of the total 
energy. Such little sensitivity to the lateral position of 
iron allows us to consider the surface essentially flat. 



Number of Fe layers 
(4 Fe per layer) 


E b (meV) 


One 


-666 


Two 


-258 


Three 


-160 


Four 


-159 


Five 


-146 



TABLE I: The dependence of the binding energy on the thick- 
ness of close-packed Fe structure with the closest layer located 
in position "1" of Figure [He) . 

Figure G2Ja) shows the binding (Ei,) and deformation 
(Edeform) energies as a function of the distance z from 
the surface. Since the AI2O3 slab experiences surface re- 
arrangement, z must be defined with respect to a fixed 
reference, in our case, the bottom of the unit cell (see 
Figure [2]). For convenience, we use a constant shift z re f 
to have the minimum of the interaction energy at a rea- 
sonable distance of 2.25 A. The strong contribution of 



energy as 



Edeform — ^ \Edeformed[Al20 3 \ 



Edeform to Ej,, shown in Figure (3Ja), is caused by the 
considerable rearrangement of the surface atoms to ac- 
commodate the absorbed iron. The origin of the non- 
monotonic variation of the deformation energy is clari- 
fied in Fig. \S[.b), in which we show the vertical shifts 
of the outer Al and O layers, z Al and z° , with respect 
to the O position when Fe is not present (z — * 00). The 
outer positively charged Al layer first rises above the sub- 
strate to be closer to the approaching Fe film and accept 
E r eiaxed[Ai 2 o 3 ]\ i (2) more charge (as has been shown previously for the Zr- 
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terminated Ni/ZrC>2 system [6 if. the charge transfer be- 
tween the outer Zr and Ni, as well as the accumulation 
of charge in the gap between them, is significant). The 
maximum rise for the Al layer happens at the interlayer 
Al-Fe distance of 2.7 A(z=2.8 A in Fig. [^b)), which 
corresponds to the kink in the deformation energy curve. 
After that point, both Al and O layers are pushed down- 
wards. For z < 1.5 A, one can also expect noticeable 
deformation of the Fe film, and the decomposition of the 
total Fe binding into the independent Fe-Fe and Fe-A^Oa 
contributions may not be accurate anymore. 
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FIG. 3: (color online), (a) The binding energy of Fe-Al2C>3 
(ab initio values with fitted Morse potential) and deformation 
energy of the AI2O3 slab as functions of the distance from the 
surface, (b) Vertical coordinates of the outer Al and O layers 
of AI2O3 as function of the distance from the surface. 



The ab initio binding energy can be decomposed as a 
sum of contributions from each atom in the Fe "column" 
(assuming that the many-body Fe-Fe effects are appro- 
priately described by an external interaction): 



Fitting gives D — 153 meV, a = 1.268 A -1 , and z = 
2.219 A. Fi gure [3ja) shows the fit of Ef, calculated as a 
sum of three Morse interactions. 

The strength of the Fe-A^Oa interaction (153 meV) 
amounts to only ^10% of the Fe-Fe interlayer binding 
energy 24-1. 42 eV, depending on the thickness of the 
Fe film), which implies that the latter is not significantly 
affected by the presence of the substrate. Hence, the as- 
sumption of decomposing the total Fe binding into the 
separate Fe-Fe and Fe-A^Oa interactions is justified a 
posteriori. Our simple Fe-A^Oa Morse interaction nat- 
urally incorporates the whole alumina surface deforma- 
tion, since the binding energy has been defined with re- 
spect to an ideal AI2O3 slab and an isolated Fe film, 
as described in equation |T]) and Figure O Within this 
framework, MD simulations on supported particles can 
be efficiently performed with reasonable accuracy by con- 
sidering the z component of the particle interaction with 
the flat substrate. 

As we are interested in the properties of binary 
Fei-zCz nanoparticles, we need to evaluate the interac- 
tion between carbon surrounded by iron and the alumina 
substrate (Fe-embedded-C-A^Oa interaction). Starting 
from three layers of adsorbed Fe, we substitute one Fe 
with one C atom bringing the concentration to x ~8.3%. 
Depending on the position in which the substitution takes 
place, the binding energy ^[FenCi-A^Oa] fluctuates 
while remaining similar to E^Fe^-A^Oa], which sug- 
gests some importance of many-body effects in the Fe- 
embedded-C-Al 2 03 interaction. Because the substrate 
represents only a small perturbation to the total bind- 
ing of the nanoparticle and the concentration of C in our 
simulations is small, we choose to treat the C-AI2O3 and 
Fe-A^Oa interactions in the same way, i.e. we use the Fe- 
AI2O3 Morse potential for both atom types. The validity 
of this approximation has been addressed by performing 
MD simulations, similar to those described in Section liTTl 
but here we vary the Morse parameters for the C-AI2O3 
interaction. Our tests for supported Fe3ooC3o nanoparti- 
cles reveal that the thermodynamics of the systems is not 
very sensitive to the particular value of the C-AI2O3 bind- 
ing. For example, increasing or decreasing the strength 
of the C-AI2O3 interaction for the Fe3ooC3o nanoparti- 
cles by a factor of two changes the melting temperature 
by 1-2 %, the amount comparable to the statistical error 
of the melting temperature determination in our simula- 
tions (see Section III). We conclude that the simple ap- 
proximation of employing the same description for the Fe 
and C interactions with the substrate is suitable for this 
study. 



(3) 



III. MOLECULAR DYNAMICS 



where the ener gy p er atom Ei is taken to be a simple 
Morse potential 62] as: 

Ei = D {exp [-2a { Zi - z )} - 2exp [-a {z t - z )]} (4) 



Among the taxonomy of thermodynamics phenomena 
for nanoparticles, melting has been subject of consider- 
able interest. The characteristics of the melting process 
depend on a variety of parameters such as size [26l . [63j 
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and shape of the particles [64| , concentration of impuri- 
ties [65[ , and presence of substrates [H, Q . 

In this section we are interested in the effect of size, 
carbon presence and AI2O3 substrate on the melting and 
solidification of binary Fe-C nanoparticles. We address 
these tasks by analyzing, with classical MD simulations, 
the liquidus and solidus lines of their phase diagrams, in 
function of the aforementioned parameters. 



A. Methods 

Molecular Dynamics Simulations. MD simulations 
are carried in the NVT ensemble using the Verlet al- 
gorithm [6(1 [6?| with a time step At = 1.0 fs. Of the 
several methods developed for controlling the temper- 
ature in MD simulations [H, IM H [Dl Izl, Berend- 
sen and Nose-Hoover thermostats are most commonly 
used. In Fig. 3] we compare the temperature disper- 
sion AT of the two thermostats with respect to that of 
the canonical distribution AT canon i ca i [73j at T = 400 
K for small nanoparticles (Np e = 50). We observe that 
the widely used Berendsen thermostat is more sensitive 
to the choice of the coupling constant. Inset in Figure 
[4] shows the distribution of instantaneous kinetic tem- 
perature for Berendsen thermostat for the typical value 
of At/r = 0.1 [74], the Nose-Hoover thermostat for r 
= 25 fs, and the canonical distribution at T = 400 K. 
The Nose-Hoover thermostat (<tt ~ 47 K) reproduces 
the canonical distribution^^ ~ 46 K) much better than 
the Berendsen thermostat (o~t ~ 17 K), making it a bet- 
ter choice for our constant temperature simulations. 

Interatomic interaction. Fe-Fe, Fe-C. and C-C in- 
teractions are described by Born-Mayer [59|, [6(3], John- 
son [75L [76| and Lennard- Jones [73 potentials, respec- 
tively. These interaction models are discussed in detail 
elsewhere [76]. The Morse potential, introduced in Sec- 
tion HH is used to model the Fe-A^Oa interaction; the 
Fe-embedded-C-Al 2 03 interaction (C is diffused in Fe) is 
taken to be identical to Fe-Al 2 03 as discussed before. 

Initial configurations. To avoid excessive temperature 
fluctuations in the MD simulations of the nanoparticles 
one should start from the most stable configurations. We 
search for the best possible energy minima by randomly 
arranging atoms in a spherical nanoparticle, carefully op- 
timizing the positions of iron and carbon atoms and fi- 
nally annealing the nanoparticles for 6 x 10 6 MD itera- 
tions (6 ns). The annealing is performed in the following 
way: the nanoparticle is first heated to high temperature 
(from 1000 K to 1400 K depending on the size of the par- 
ticle) for 0.6 x 10 6 steps, kept at constant temperature 
for another 0.6 x 10 6 iterations, and finally cooled to K 
during the remaining 4.8 x 10 6 MD steps. 

Definition and determination of melting temperature. 
In our work, the melting phenomenon is analyzed by per- 
forming several MD simulations starting at about 300 K 
below the expected melting point with temperature incre- 
ments of 10 K for small [N < 100) and 20 K for large clus- 
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FIG. 4: (color online). Comparison of Berendsen and Nose- 
Hoover thermostat for Nfc = 50 at 400 K for At — 1 fs: tem- 
perature dispersion for Nose-Hoover is closer to the canonical 
dispersion. Inset shows the distribution of instantaneous ki- 
netic temperature for Berendsen thermostat for At/r — 0.1, 
Nose-Hoover thermostat for r = 25 fs, and the canonical dis- 
tribution. 



ters (with 5 K upon approaching the transition). Only 
the lowest temperature simulations begin from the an- 
nealed initial structures: the others start from the final 
configurations (positions, forces, velocities) of the pre- 
ceding temperature simulation. Data gathering of the 
energies and other averages are performed over 10 6 MD 
steps. 

Several dynamical and structural properties such as 
total energy, Lindcmann index, diffusion coefficients, and 
pair correlation functions can be used to identify phase 
transitions in nanoparticles [zl, [79|. Here, melting is 
characterized by the temperature dependence of the total 
energy (caloric curve), by the change in the total energy 
with time, and by the variation of the Lindemann index 
with respect to temperature 80]. The Lindemann index 
S represents the root-mean-square relative bond-length 
fluctuation: 



f V 

N(N - 1) ^ 

v ' »<j 



(5) 



where ry is the distance between atom i and j, N is the 
number of particles and the average is calculated over 
an MD run at a given T. The melting point, which 
defines the temperature at which a solid becomes liq- 
uid, is a macroscopic concept for pure and bulk sys- 
tems. Both finite-size and presence of more than one 
atomic species make the melting transition a continuous 
phenomenon that occurs over a range of temperatures, 
AT m , in which solid and liquid phases coexist with dif- 
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FIG. 5: Estimation of the melting temperature T m as max- 
solid point for Fe4oo nanoparticle from the caloric curve (panel 
a) and from the Lindemann index plot (panel b). 



ferent fractions [8l|, |82j, |83J] . To have a specific value of 
T m instead of a range, we define the melting tempera- 
ture T m as the "max-solid poinf which represents the 
maximum temperature at which the solid and the liquid 
phases coexist (the locus of all the max-solid points is 
the liquidus). Above T m , no solid phase is present. Note 
that within this definition of T m , we also identify plastic- 
viscous nanoparticles as "liquid" [84J . Figure \5\ shows an 
example of T m calculated from the caloric curve (panel a) 
and from the Lindemann index (panel b). Similarly, the 
min-liquid point is the minimum temperature at which 
the solid and the liquid phases coexist (the locus of all 
the min-liquid points is the solidus). The difference be- 
tween the energies of the particle at the max-solid and 
at the min-liquid points defines the enthalpy of melting 
AH m . 



B. Melting of nanoparticles 

With the aforementioned method, we investigate pure 
nanoparticles of size Np e = 80 — 1000 atoms (diam- 
eter d ~ 1-3 nm). Caloric curves for particles with 
Np e > 100 show small melting intervals (AT m < 10 K). 
For smaller clusters (Np e ~ 80 — 100 atoms, d ~1 nm) 
the characterization becomes difficult because both the 
caloric curves and the Lindemann indices fluctuate over 
wide intervals of temperatures in which the liquid and 
solid phases coexist in dynamic equilibrium. This well 
known phenomenon, called dynamic coexistence melting 
[78| and shown in Figures [U^a) and Efb), is caused by 
the multitude of different metastable solid phases present 
at the nanoscale, which have similar free energies, simi- 
lar volumes (at constant pressure volumes are allowed to 
change) and different surface arrangements. The parti- 
cles are quasi-plastic by continuously changing their state 
while alternating metastable configurations and liquid 
states (bi-stability). Thus, the observations of thermal 
properties inside the dynamic coexistence melting inter- 
val AT m are affected by two types of fluctuations: physi- 
cal, due to coexistence of phases (cannot be avoided) , and 
statistical (can be reduced by increasing the total time 
of the MD simulation). In particular, the fluctuations of 
the Lindemann index can be captured by analyzing the 
standard deviation o~li. Figure [He) illustrates the phe- 
nomenon: by plotting o~n versus T, we can estimate the 
dynamic coexistence melting interval AT m , the max-solid 
and min-liquid points. Pure liquid and solid temperature 
ranges will be the ones with negligible <tli- 

Figure [7J shows T m of pure Fe nanoparticles in the 
whole range of sizes (Np e = 80 — 1000) for free and sup- 
ported clusters. Our results, in agreement with other 
theoretical [H, [33. computational (l3l . [3^ . and exper- 
imental studies [20, H3, HI, [H, H3|, predict a decrease 
in the melting temperature inversely proportional to the 
cluster diameter [85|]. The behavior of T m can be de- 
scribed by the model based on the Gibbs-Thomson equa- 
tion [2H HH, H|| in function of bulk melting temperature 
T^ lk , effective diameter of the particle d, latent heat of 
melting AH SV , and solid- vapor interfacial energy 7 S „ (87j • 

The melting point of bulk Fe, obtained by extrapolat- 
ing the fit to d -> 00 (T^(N Fe = 00) -1416 K) is 
— 20% below the real one [881 ]. indicating that our Fe-Fe 
Born Mayer many body potential is slightly undcrbind- 
ing. The systematic shift in melting temperatures should 
not affect the main conclusions of our study since we are 
interested in the trends of the liquidus lines in the phase 
diagrams rather than their precise values. The supported 
particles considered here have higher melting point than 
that of the free clusters. In fact, the attractive interac- 
tion with the substrate (AI2O3) induces flattening of the 
particle and increases the effective diameter in agreement 
with previous studies [l5| . 

Melting temperatures for magic-size clusters (N = 
55,147,309,561,923) in our range of simulations are 
shown in the inset of Figure [7] Due to their inherent 
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FIG. 6: Melting phenomenon for small nanoparticles 
(Nfe ~80-100 atoms, d ~1 nm): a) caloric curve, b) Lin- 
demann index with respect to temperature, c) standard devi- 
ation a li of the Lindemann index to identify the max-solid 
and min-liquid points. 



1400 
1300 
1200 
1100 
1000 
900 
800 





1300 


■ s923 








1200 

„ 1100 
& E 1000 


■ s 5v561 

^*s309 

"Vs. 147 






" 1000'~v^< 

800 *\\ 
600 


H 900 
800 








700 










0.3 0.4 0.5 0.6 0.7 0.8 


0.8 1.0 


400 




d" 1 (nrn" 1 ) 








• 










200 155"-X 










100 *-.«\ 









i 


80 V " 



i 




0.3 0.4 0.5 


0.6 


0.7 0.8 0.9 


1.0 





d" 1 (nm" 1 ) 

FIG. 7: (color online). Melting temperature versus the inverse 
of particle diameters. Free non-magic sizes (•) and their lin- 
ear fit (dashed line); supported non-magic sizes (A) and their 
linear fit (solid line); Supported clusters have higher melting 
point due to the decreased curvature of the surface [TBI] . Melt- 
ing temperatures of free magic sizes (o) and supported magic 
sizes (A) clusters were compared with the linear fit lines of 
non-magic sizes (inset). 



magic-size if free, is not magic if supported (because of 
the shape and sometimes structural modifications due to 
interaction with the attractive surface; the phenomenon 
is addressed in Section pV[) ) . In summary, the small sup- 
ported magic-size clusters have lower melting point than 
the corresponding unsupported ones. 

The melting temperature of very small clusters (N < 
80) is a function not only of their size but also of their spe- 
cific structure. It has been shown theoretically that melt- 
ing temperatures of clusters with several or tens of atoms 
can be abnormally high, e ven above the corresponding 
bulk values [§(| HH IS ■ This phenomenon is related 
to atomistic processes of structure isomerization (geom- 
etry reconstruction) [8^, [9(| HH [H[ or electronic struc- 
ture change (formation of strong covalent bonds) [92], [94[ . 
Since the thermodynamics of magic-size nanoparticles is 
significantly influenced by their structure, they experi- 
ence unusual peculiarities. In this Section we focus our 
analysis only on the non-magic size ones. 



Phase diagram of free and supported Fe-C 
nanoparticles 



symmetry (icosahedral or decahedral for small clusters in 
our case), free magic-size clusters (smaller than Np e = 
309) are very stable with melting temperatures higher 
than that of the non-magic ones. This does not apply 
for supported clusters. In fact, the presence of substrate 
interaction changes the magic-size sequence: a cluster of 



Analysis of phase diagrams. To understand how in- 
clusion of carbon atoms influences the thermal behavior 
of the catalyst nanoparticles, we determine the melting 
temperatures as a function of carbon concentration x 
ranging from zero to up to <~ 16 %. To appropriately 
model the nanocatalyst in a nanotube growth process, 
instead of substituting Fe with C to increase the concen- 
tration of carbon, we add Nc atoms to the Fe particle. 
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Since Np e remains constant and Nc increases, the con- 
centration of carbon is defined as x c = Nc/(Np e + Nc)- 
We plot the locus of the max-solid points which repre- 
sent the liquidus of the Fe-C phase diagram. In small 
Fe nanoparticles one might expect the average radius to 
increase noticeably with the addition of a few C atoms 
(r 3 oc N) and hence their melting temperature would 



also increase ((Tj 



bulk 



T m )/T; 



bulk 



d 1 ). However, our 



tests show that addition of C (x c < 16%) does not signif- 
icantly change the volume of the nanoparticle indicating 
that C behaves as an interstitial solute in Fe nanoparti- 
cles as it does in bulk Fe (95j . 

Figure [5] shows phase diagrams (T m versus x c in 
atomic %) for free and Al 2 03-supported particles with 
Np e = 80, 100, and 200 based on caloric curve and Lin- 
demann index analysis. All the data sets show a similar 
trend in function of C concentration: T m decreases al- 
most linearly at low x c and then increases for all the 
higher x c considered. The exact functional form is diffi- 
cult to determine because of the dispersion in the data, 
however the observed "V'-shape dependence is consistent 
with that in the bulk Fe-C phase diagram [8q . Hence, 
by using the least square method we approximate the liq- 
uidus with a set of two straight lines, the intersection of 
which gives the eutectic point (x^ ut , T eut ) [96]. This pro- 
cedure allows us to estimate this invariant point with an 
accuracy of 1 % and 12 K for x^ ut and T eut , respectively. 

We observe that as the particle size is reduced, the eu- 
tectic point for free and supported nanoparticles moves 
toward lower temperatures and lower concentrations, in- 
dicating that the solubility of C decreases as well. To the 
best of our knowledge this phenomenon has never been 
reported before. Explanation of its origin will require 
a more detailed study of the behavior of dissolved C at 
various concentrations (changes in the distribution of C 
across the particle, possible formation of stable carbides 
etc.). 

Implications for carbon nanotube growth. The "V'- 
shape liquidus feature of Fe-C nanoparticles observed in 
our simulations allows for the VLS interpretation of ex- 
perimental results for catalytic activity of Fe. For a given 
temperature above the eutectic temperature T eut , the 
dissolution of carbon in a metal catalyst initially induces 
the liquefaction of the particle by lowering its melting 
temperature [111, H3] • As the catalyst becomes less vis- 
cous, the diffusion of carbon in the particle increases. 
High catalytic activity has been observed with associ- 
ated liquefied particles during the growth of SWNT by 
CVD method [97| . Some experiments have shown pres- 
ence of liquid-like features and liquid layers on nanopar- 
ticles [H, |99[ before complete melting. These features 
would enhance the diffusion of carbon and the subsequent 
melting of the particles. In our case, the nanoparticles are 
so small that the surface effects on the melting tempera- 
ture should be dominant. In fact, one cannot distinguish 
between surface and bulk layers in a particle with 200 
atoms of radius ~ 0.8 nm, which is approximately three 
close-packed layers. 
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FIG. 8: (color online). Fe-C phase diagrams obtained with ad- 
dition of C (upto ~ 16%) to particles with iV^e =80,100, and 
200: a) free nanoparticles and b) AI2O3 supported nanopar- 
ticles. 



According to the data in Fig. [8j if during CVD ex- 
periments the size of the catalyst particle were in a given 
range (d ~ 1 — 2 nm), the smaller Fe nanoparticles (d ~ 1 
nm) would liquefy earlier then the bigger ones (d ~ 2 
nm) . Consequently, if liquefaction was a prerequisite for 
catalytic activity, the small particles would begin to pro- 
duce nanotubes earlier. Experimental verification of this 
hypothesis is challenging: very small catalyst nanoparti- 
cles (d =1 nm) have a fast rate of coalescence durin g the 



looLfToil 



reaction and coagulate to form bigger clusters 
Following the VLS model, the growth of very small nan- 
otubes using metallic catalyst would be possible only if 
the reaction temperature were above T eut (to liquefy the 
particle) and below the temperature at which particles 
begin to coalesce. Hence, controlling the diameter of nan- 
otubes grown with CVD in the small size range might be 
difficult. So far, the smallest reported nanotube (d = 4 
A, the most internal tube in a multiwalled nanotube) has 
been grown by a rc-discharge method without any metal- 
lic catalyst fl02] . 
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In the x c > Xg Ut region, our results also show that once 
the dissolved carbon concentration reaches a point where 
the corresponding melting temperature (from the phase 
diagram) exceeds the reaction temperature, the parti- 
cle starts solidifying as a two-phase system composed of 
solid carbide and Fe-rich liquid. This solidification grad- 
ually reduces the average surface mobility of the catalytic 
species (Fe) and might affect the catalytic activity of the 
nanoparticle. During the process of nanotube growth, 
conditions such as sufficient carbon concentration and 
temperature can induce the formation of stable carbides. 
Cementite (FeaC) and other iron carbides have been ob- 
served in experiments after the reactio n of hydrocar bon 
or carbon monoxide with iron catalysts [l03lll04ill05j |. It 
has been reported that such stable carbides act a s poison 
by terminating the growth of nanotube [HI, Il06l | . 



IV. 



PECULIARITIES OF SMALL FREE AND 
SUPPORTED FE CLUSTERS 



In this section we analyze the structural properties of 
pure small Fe clusters, free and supported on AI2O3. We 
find that as the size decreas es Fe clus ters assume sta- 
ble polyhedron configurations [l07i Il08j |. The same MD 
technique described in the previous section is adapted 
for searching the stable configurations of small clusters 
(N ~ 50 — 80). Twenty random configurations are gen- 
erated per each size and optimized by annealing. The 
whole annealing process contains 10 7 MD iterations (10 
ns). Free and supported clusters are heated to 1200 K 
and 1400 K, respectively, before being slowly cooled to 
K. Among the twenty annealed structures, we consider 
the lowest energy configuration as the global minimum 
for each size. For free clusters with N < 72, the lowest 
5 energies are nearly degenerate (inset in Fig. [5]) which 
validates our approach. For bigger particles the found 
minima might be sub-optimal as the difficulty of finding 
the global minima rapidly increases. 

The minimum energies per atom are fitted to the 
Eave = a + bN^ 1 / 3 dependence, where, in the case of 
spherical particles, the parameters a and b would repre- 
sent the bulk energy per atom and the surface-cr e atio n 
destabilization energy cost, respectively [l09l ll IOL llllj . 
Our fit uses only two parameters (a, b) instead of the 
four of the mor e general form ula E ave — a + 6A~ 1 / 3 + 
cN-W + dN' 1 [lOl, [HO, [HU to avoid overfitting due to 
the limited number of energy points. We include the best 
available minima obtained with the aforementioned pro- 
cedure for bigger (N = 100, 120, 160 and 200) clusters to 
ensure the correct assymptotic behavior for the nanopar- 
ticles of large sizes. To capture the size dependence only 
we exlude sizes N ~ 55 (N — 54, 55, and 56) from the 
fit, because these nanoparticls acheive additional stabil- 
ity by forming highly symmetric icosahedra (see Fig. [5]). 
We obtain a free = -4.29 ± 0.01 eV, b free = 2.42 ± 0.03 



values of a are very close to the bulk fee cohesive en- 
ergy of -4.29 eV/atom; the two values of b slightly differ 
because of the reduced total energy cost to create and 
modify surface for the supported cluster [l5| . 

Figures [10] and QT] illustrate the differences between 
the calculated and the fitted energies for free and sup- 
ported clusters. The prominent negative peaks in Fig. 
[TO] for sizes N = 55,61,64,71,75, and 77 correspond to 
the magic size structures in our sequence. The evolution 
of the nanoparticle configuration with size is depicted in 
the same picture. By comparing Figures [TO] with [11] we 
observe that the attractive substrate has little or no effect 
on the internal structure of small clusters (N < 70) . Due 
to the spherical (or nearly-spherical) arrangements, such 
small icosahedra are not significantly deformed by attrac- 
tive substrates, unless the adsorption potentials is com- 
parable to the internal atomic binding energy of the clus- 
ter. Big supported particles behaves differently. In fact, 
as the size increases beyond N > 70, alumina-supported 
clusters form fee arrangements (Fig. [TTj) with the {111} 
planes at the particle/substrate interface. Elsewhere the 
particles have simple close-packed facets or local arrange- 
ments of them. 

Big unsupported particles are different than the sup- 
ported ones. Unsupported clusters of sizes between N ~ 
70 and N ~ 200 have decahedra or icosah edra arrange- 
ments (multi- twinned structures [ll2[ Ill3| h Decahedra 
are more frequent for medium size particles (N ~70 to 
~100) while icosahedra appear more often for bigger sizes 
(N ^100 to ^200). Figure [TO] shows three particular ex- 
amples of decahedra clusters with TV = 71, 75, and 77. 
As the size exceeds N ~ 200, the particles approach the 
bulk configuration and the fee and icosahedra structures 
tend to become degenerate. For instance, the cluster at 
N = 500 is a good example of an fee particle with ap- 
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FIG. 9: (color online). The minimum energies of free clusters 
and the fitting curve using E ave = a + 6JV -1 / 3 . The 5 lowest 
energies are shown in the insert for cluster's sizes from 50 to 
80. 
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FIG. 10: Energies of best minima for free clusters relative 
to Eave, with a selection of stable structures. For free clus- 
ters, the polyhedron^fcc transition is continuous and occurs 
around N ~ 200 - 500 (not shown). 



clusters. However, for supported particles, the fee config- 
uration with flat facets (along the directions of minimum 
surface energy given by the Wulff plots [1171 [ll8() has 
a larger contact area with the substrate than the poly- 
hedra do. Therefore, the attractive substrate interaction 
reduces the surface/interface energy of fee the most, re- 
sulting in a lower critical r (even though the overall sur- 
face area of the fee structures could be larger than that 
of the icosahedra). 

Our simulations reveal that the substrate attractive in- 
teraction tends to preserve the bulk structure for small 
particles. This demonstrates that although the Fe-AlsOs 
interaction is relatively weak, it is able to influence the 
balance between the bulk and surface energies in compet- 
ing configurations and utlimately determine the particle's 
ground state. One can also expect to observe this effect 
in simulations with different models of the Fe-Fe and Fe- 
AI2O3 interactions, as long as their relative strength is 
comparable with ours. 
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FIG. 11: Energies of best minima for supported clusters rela- 
tive to Eave, with a selection of structures. For alumina sup- 
ported clusters, the polyhedron^fcc transition occurs around 
TV ~ 70. 



propriate close-packed faceting and minimal surface-edge 
reorganization. 

The polyhedron— >fcc transition that occurs at lower 
sizes for supported clusters can be explaine d in term s 
of surface/interface and elastic strain energies 
In a groundstate configuration, the former (latter) en- 
ergy dominates over the other for a small (big) particle. 
Particles with radii of r < rp (m is the critical size for 
the polyhedron-fee transition [116j |) tend to form polyhe- 
dral structures (icosahedra or decahedra) to minimize the 
surface/interface energy, while clusters with r > ro will 
prefer the bulk configuration (fee) to minimize the elas- 
tic strain energy |ll5| . This is true for free and supported 



V. CONCLUSIONS 

In this paper we have investigated the behavior of free 
and alumina-supported Fe-C nanoparticles. We observe 
interesting phenomena that can be attributed to the pres- 
ence of the substrate. The main results of the present 
study can be summarized as follows: 

(i) The total Fe binding can be conveniently decom- 
posed into independent Fe-Fe and Fe-A^Os parts: ac- 
cording to our ab initio calculations the two differ by 
about an order of magnitude. Moreover, the corrugation 
of the Fe-Al203 interaction is much smaller than the av- 
erage adsorption energy. This allows us to parameterize 
the Fe-A^Oa interaction as a Morse potential, which in- 
cludes the deformation energy of the substrate. 

(ii) The thermal behavior of pure-Fe particles is sim- 
ulated with classical MD techniques. We observe the 
reduction of melting temperature as a function of the di- 
ameter, in agreement with the Gibbs-Thomson law. We 
also show that supported particles have higher melting 
points than the unsupported ones. 

(iii) We calculate the liquidi on the phase diagrams for 
a range of Fe-C nanoparticles (d ~ 1.1 — 1.6 nm) and show 
that they are characterized by the presence of eutectic 
points in which the eutectic concentration depends on 
the size of the particles. These phenomena may have 
important effects on the growth of carbon nanotubes by 
CVD, as discussed in Section III C. 

(iv) We find that the optimized configurations of very 
small pure Fe clusters (down to Np e ~ 50) have icosa- 
hedron or decahedron structures. Bigger clusters tend to 
have close-packed configurations. We show that the size 
at which the cluster undergoes the polyhedron— >close- 
packed transition depends on the substrate. In particu- 
lar, AI2O3 supported Fe clusters have n£.' 2 ° 3 - 70 while 



free cluster have N; 



free 



200 - 500. 
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The thermodynamic study of nanoparticles will be ex- 
tended with future ab initio characterizations carried out 
in this laboratory, and with experimental investigations 
performed by our collaborators. This work stimulates 
more comprehensive studies of the role of substrates on 
the thermodynamics of small particles, to understand the 
fundamental factors controlling the catalytic properties 



of small clusters. 
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